DATA_SECTION
  init_int styr
  init_int endyr

  init_int nobs_F
  init_matrix Fdat(1,nobs_F,1,4)

  init_int nobs_S
  init_matrix Sdat(1,nobs_S,1,4)

  vector tend_F(1,nobs_F)
  vector age_mid_F(1,nobs_F)
  vector cat_F(1,nobs_F)
  vector coh_F(1,nobs_F)
  ivector yr_F(1,nobs_F)
  ivector mo_F(1,nobs_F)

  vector tend_S(1,nobs_S)
  vector age_mid_S(1,nobs_S)
  vector cat_S(1,nobs_S)
  vector coh_S(1,nobs_S)
  ivector yr_S(1,nobs_S)
  ivector mo_S(1,nobs_S)

  int minK
  int maxK
  int ncoh
  int nages
  number offset
  number M
 LOCAL_CALCS
  // Get data into vectors (from matrix
  tend_F    = column(Fdat,1);
  age_mid_F = column(Fdat,2);
  cat_F     = column(Fdat,3);
  coh_F     = column(Fdat,4);
  for (int i=1;i<=nobs_F;i++)
  {
    yr_F(i) = int(tend_F(i));    
    mo_F(i) = 12*(tend_F(i) - yr_F(i)+.01);    
  }

  tend_S    = column(Sdat,1);
  age_mid_S = column(Sdat,2);
  cat_S     = column(Sdat,3);
  coh_S     = column(Sdat,4);
  for (int i=1;i<=nobs_S;i++)
  {
    yr_S(i) = int(tend_S(i));    
    mo_S(i) = 12*(tend_S(i) - yr_S(i)+.01);    
  }
  
  // Compute the number of observations by 
  maxK = max(max(coh_S),max(coh_F));
  minK = min(min(coh_S),min(coh_F));
  ncoh = maxK-minK+1;
  nages = 10;
  M=0.2;
   offset = 0.0;
 END_CALCS

INITIALIZATION_SECTION
  log_N0 8.;
  log_F -2.;

PARAMETER_SECTION
  init_vector log_N0(1,ncoh);
  init_number log_F
  // matrix S_trans(1,ntyp+1,1,ntyp+1) // Survival transition matrix from M and each fishery (or survey)
  vector N0(1,ncoh);
  matrix prC(minK,maxK,1,nages) 
  matrix prO(minK,maxK,1,nages) 
  matrix resid(minK,maxK,1,nages)

  // By cohort     
  matrix Nc(1,ncoh,1,nages)    
  matrix Sc(1,ncoh,1,nages)    
  matrix Zc(1,ncoh,1,nages)    
  matrix Fc(1,ncoh,1,nages)    
  matrix Cc(1,ncoh,1,nages)    

  // Flat structure
  matrix N(styr,endyr,1,nages)    
  matrix S(styr,endyr,1,nages)    
  matrix F(styr,endyr,1,nages)    
  matrix Z(styr,endyr,1,nages)    
  matrix C(styr,endyr,1,nages)    

  number notR
  number totR
  number pr_notR
  number totprR  
  objective_function_value nll;

PROCEDURE_SECTION
  N0 = mfexp(log_N0);
  F = mfexp(log_F);
  Z = F + M;
  S = mfexp(-Z);
  for (int i=styr; i<=endyr; i++)
    N(i,1) = N0(i);
  for (int j=2;j<=nages;j++)
    N(1,j) = N(1,j-1)*S(1,j-1);
  N(1,nages) = N(1,nages)/(1.-S(1,nages));
  for (int i=2;i<=ncoh;i++)
  {
    N(i)(2,nages) = ++elem_prod(S(i-1),N(i-1))(1,nages-1);
    N(i,nages)   +=  N(i-1,nages)*S(i-1,nages);
  }
  C = elem_prod(elem_div(elem_prod(N,F),Z),1.-S);

  nll = 12.5*norm2(F-.1) + 50.*norm2(C-20.) + 50.*norm2(log_N0-mean(log_N0));

REPORT_SECTION
  report <<"N"<<endl;
  report << N <<endl;
  report <<"F"<<endl;
  report << F <<endl;
  report <<"C"<<endl;
  report << C <<endl;
